I 


AD-A096  548  DEFENCE  RESEARCH  ESTABLISHMENT  ATLANTIC  DARTMOUTH  (NO— ETC  F/6  20/1 

N0RM2L :  AN  INTERACTIVE  COMPUTER  PROGRAM  FOR  ACOUSTIC  NORMAL  MOD~ETC<U) 
DEC  80  D  D  ELLIS 

UNCLASSIFIED  0REA-TM-80/K  NL 


FILE  COPY  AD  AO 96 54  8 


NORM2L 

AN  INTERACTIVE  COMPUTER  PROGRAM 
FOR  ACOUSTIC  NORMAL  MODE  CALCULATIONS 
FOR  THE  PEKERIS  MODEL 


DEFENCE  RESEARCH  ESTABLISHMENT 
ATLANTIC 

DARTMOUTH  N.S. 


O.R.E.A.  TECHNICAL  MEMORANDUM  80/K 

f:  ' 


'  t  ‘  '  , 

\ 

~~ . 

J?-  NORM2L 

_  -  AN  INTERACTIVE  COMPUTER  PROGRAM 
FOR  ACOUSTIC  NORMAL  MODE  CALCULATIONS 
FOR  THE  PEKERIS  MODEL, 


O.  D.  kllis 


80  i 


Approved  by  R.  F  BrOWfl  Director  /  Under woter  Acoustics  division 


DISTRIBUTION  approved  by 


CJ 


CHIEF  D  R  E  A 


RESEARCH  AND  DEVELOPMENT  BRANCH 

DEPARTMENT  OF  NATIONAL  DEFENCE 

CANADA 


ABSTRACT 


The  interactive  computer  program,  N0RM2L,  calculates  the 
discrete  normal  modes  and  acoustic  propagation  loss  for  the  Pekeris 
model  of  the  ocean.  The  Pekeris  model  is  a  simple  two-layer  model 
in  which  the  two  layers  represent  the  seawater  and  seabed.  For  many 
shallow-water  environments,  the  model  is  a  reasonable  approximation  to 
the  actual  physical  situation  and  can  be  used  to  investigate  acoustic 
propagation  at  low  frequencies. 

For  ease  of  future  expansion  and  modification,  the  program 
N0RM2L  is  written  in  modular  form  in  FORTRAN.  The  results  of  N0RM2L 
are  compared  with  those  of  other  computer  programs. 
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RESUME 


Le  programme  d'ordinateur  interactif  N0RM2L  sert  A  calculer  les 
modes  normaux  discrets  et  la  perte  de  propagation  acoustique  pour  le 
modele  de  Pekeris  de  I'ocdan.  Le  moddle  de  Pekeris  est  un  module  simple 
A  deux  couches  dans  lequel  les  deux  couches  reprdsentent  l’eau  de  mer  et 
le  fond  de  la  mer.  Pour  plusieurs  environnements  en  eau  peu  profonde, 
ce  mod&le  est  une  approximation  raisonnable  de  la  situation  rdelle  et 
peut  servir  5  dtudier  la  propagation  des  ondes  acoustiques  aux  basses 
frequences . 

Pour  en  faciliter  l'expansion  et  la  modification  futures,  le 
programme  N0RM2L  est  dcrit  sous  forme  modulaire  en  langage  FORTRAN. 

L' auteur  compare  les  rdsultats  du  N0RM2L  A  ceux  des  autres  programmes 
d' ordinateur. 
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INTRODUCTION 


NORM2L  is  an  interactive  computer  program  which 
calculates  the  discrete  normal  modes  and  propagation  loss 
for  the  Pekeris  [1]  acoustic  model  described  in  section  2.1. 
Although  this  model  is  simple,  it  contains  all  the  essential 
features  of  normal  mode  theory,  and  for  many  shallow  water 
environments  is  a  reasonable  approximation  to  the  actual 
physical  situation. 

The  main  reason  for  writing  N0RM2L  was  to  have  a 
convenient  program  for  becoming  acquainted  with  normal  mode 
calculations.  The  interactive  program  and  simple  model  allow 
rough  answers  to  be  obtained  quickly  and  cheaply.  Since 
many  of  the  solutions  are  known  analytically,  the  extensive 
calculations  of  a  general  normal  mode  program  are  not 
necessary.  Also,  the  interactive  nature  of  the  program 
allows  the  user  to  see  the  results  immediately.  The  user 
can  then  change  the  input  parameters  and  look  at  the  new  set 
of  results.  In  this  way,  the  user  can  quickly  get  a  feeling 
for  the  importance  of  the  physical  parameters. 

Another  reason  for  writing  the  program  was  to  test 
new  ideas  and  algorithms  in  normal  mode  theory  in  an  easy  to 
modify  package.  For  example,  to  help  verify  that  a 
numerical  procedure  is  working  or  that  a  particular  normal 
mode  program  has  been  implemented  properly,  the  results  from 
the  analytic  solutions  of  the  Pekeris  model  can  be  compared 
with  similar  results  from  other  normal  mode  programs  which 


obtain  their  results  numerically.  In  one  application  a 
slightly  modified  version  of  N0RM2L  was  found  valuable  in 
testing  and  debugging  a  normal  mode  program  based  on  a 
two-ended  shooting  method  [2,31.  In  another  application, 
not  foreseen  at  the  time  N0RM2L  was  written,  the  program  was 
modified  and  some  of  the  output  used  for  mode  enhancement  in 
a  preliminary  analysis  of  some  real  data  [4].  The  present 
report  does  not  discuss  these  applications,  but  restricts 
itself  to  a  description  of  the  program. 

The  main  body  of  the  report  contains  a  general 
description  of  the  model  and  program,  with  details  relegated 
to  the  appendices.  Section  2  contains  the  theoretical  basis 
with  a  summary  of  the  Pekeris  model  and  a  description  of  the 
algorithm  used  to  calculate  the  normal  mode  wave  numbers. 
Section  3  contains  a  general  description  of  the  program,  its 
structure  and  important  variables,  and  results  of  testing 
and  verification.  A  number  of  tables  are  presented  on  wave 
numbers,  mode  functions,  and  propagation  loss  to  be  used  as 
reference  for  checking  the  accuracy  of  new  or  other  normal 
mode  programs.  In  Section  4,  the  recommendations  contain 
some  possible  future  extensions  and  uses  of  the  program. 

Appendices  A  to  D  contain  further  details  on  the 
program.  Appendix  A  is  a  user's  guide  for  N0RM2L  with 
details  on  the  hardware-software  environment,  input  and 
output,  operational  instructions,  restrictions  and 
limitations,  functions  and  subroutines,  accuracy,  memory 


requirements,  and  execution  times.  Appendices  B  and  C 
present  sample  input /output,  from  the  terminal  and  line 
printer,  and  Appendix  D  contains  a  program  listing. 

Appendices  F  and  F  are  forward-looking,  giving  some 
anticipated  modifications  and  extensions  and  a  place  for  the 
user  to  add  his  comments  and  suggestions.  Some  of  these 
changes  will  be  incorporated  in  future  versions  of  N0RM2L. 
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THEORETICAL  BASIS 


2. 1  The  Pekeris  Model 

Figure  1  illustrates  the  environment  described  by 
the  Pekeris  [1]  acoustic  normal  mode  model.  The  model  has  a 
homogeneous  layer  of  thickness  h  bounded  above  by  a  pressure 
release  surface  and  bounded  below  by  a  homogeneous 
semi-infinite  half-space  of  greater  sound  speed.  The  sound 
speeds  c-|  and  Cp  and  the  densities  p  ^  and  p  p  are  constants. 
To  make  the  environment  more  realistic,  the  bottom  layer  has 
been  given  an  absorption  coefficient  a p.  A  source  at  depth 
zs  in  the  first  layer  is  operating  at  frequency  f.  The 
problem  is  to  calculate  the  propagation  loss  between  the 
source  and  a  receiver  at  depth  zf  and  horizontal  distance  r 
from  the  source. 

The  greater  sound  speed  of  the  bottom  causes  some  of 
the  acoustic  energy  to  be  totally  reflected  from  the  bottom. 
The  successive  reflections  from  the  surface  and  bottom  cause 
the  water  layer  to  act  as  a  waveguide,  and  a  standing  wave 
pattern  is  set  up  in  the  depth  direction.  These  standing 
waves  are  the  normal  modes.  Near  the  source  some  of  the 
energy  :;t.  ri  kos  the  bottom  at  angles  less  than  the  critical 
angle  and  is  partially  transmitted.  This  energy  corresponds 
to  the  continuous  spectrum  of  normal  mode  theory.  Since 
this  energy  leaks  out  to  the  bottom  after  several  bottom 
reflections,  it  is  nearly  always  neglected  in  normal  mode 
calculations.  Like  most  normal  mode  programs.  N0RM2L 
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calculates  only  the  discrete  normal  modes.  This  means  that 
the  model  is  not  applicable  when  the  source  and  receiver  are 
separated  by  horizontal  distances  less  than  several  water 
depths . 

A  standard  reference  for  underwater  acoustics  is 
Tolstoy  and  Clay  [5],  although  the  formulae  for  the  normal 
mode  equations  are  scattered  throughout  the  book.  Chapter  9 
of  the  newer  book  by  Clay  and  Medwin  [6]  contains  a  fairly 
clear  introduction  to  normal  mode  theory.  Only  the  essential 
equations  are  presented  here.  Note  that  the  formulation  is 
in  terms  of  the  velocity  potential,  and  that  only  the 
discrete  modes  are  calculated. 

For  a  frequency  f  (angular  frequency  oj  =  0^f)  the 
normal  modes  are  the  eigenfunctions  un(z)  of  the  following 
differential  equation 

d2u  2 

— -r  +  \^2  -  k2  JU  (z)  =  0  .  (  1  ) 

oz  c  n  n 

Once  the  eigenvalue  kn  is  determined  (see  below),  the 
solutions  un  for  the  two  layer  model  are  given  as  follows: 

u  (z)  -  A  sin(>,  z)  ,  o*-z<h  (2a) 

n  n  I  n  —  — 


u  (.z)  =  (p./p,)  A  sin(  y .  h)exp[-y7  (z-h)  ]  ,  z>h  (2b) 

n  1  l  n  in  zn 


where 

Yin  =  WU\  ~  kn]1/2 


( 3a ) 


(3b) 


Y-,  =  [kz  -  u‘/c 

2n  n  2 


n2il/2 


Hid 


-2  1  H  Sl"2(l,lnh). 

A„  '  2  Blh(1  -  2vlnh  n2  72„h 


(4) 


The  normalization  factor  An  is  determined  from: 

p(z)un(z)um(z)dz  =  6nm 
'  o 

where  >S  is  the  Kronecker  delta  function. 

The  conditions  of  continuity  of  pressure  and 
particle  velocity  require  that  at  the  interface  z  =  h. : 

p  u  (h  )  =  p.u  (h+) 

In  2  n 

du  (h  )  du  (h+) 

n _  _  n 

dz  dz 


(5) 


(6a) 


(6b) 


torn:  it  ion  6a  has  already  been  incorporated  into  the  equation 
for  un .  Condition  6b  gives  the  equation  used  to  determine 


tan 


f  h  [u>2 /c2 


k2|l/2,  ,  _  !2|“2_/c(  -  k»,l/2 


1  k‘ 
n 


,2/c 


(7) 


i:l  ion  of  this  equation  leads  to  a  set  of  discrete  values 


or 


u)/Ci  >  k  >  k2  >  ^  kN  >  w/c2 


(8) 


If  there  is  attenuation  in  the  bottom,  the  wave 
number  of  the  sound  speed  in  the  bottom  can  be  written  as 
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the  complex  number  kg  =  u/cp  +  ie.  Usually  e  is  quite  small 
in  comparison  to  w/cp  and  perturbation  theory  can  be  used 
[71.  Then , 


k  k  +  i  6 
n  n  n 


(9) 


wnere  , 

u'P r°° 

6  =  e  — r— -  i  I  u  ( z)  I  2  dz 

n  C2kn  Jh  n 

For  the  two  layer  model  the  integral  can  be  evaluated, 
giving 

°  A  2 

cop  ;AZ 


(10) 


r  1  n 

<5  =  e  ^ ^ - 

n  2p 0c „k  y. 


(11) 


2  2  n  2n 

Note  that  the  absorption  coefficient  ap  is  usually  given  as 
dB  per  unit  length  referenced  to  1  kHz.  Since  the 
attenuation  is  assumed  to  be  linearly  dependent  on  the 
frequency,  then  e  is  given  by 


e  =  a2  (f/1000)  ( 5 n  10)/20  (12) 

and  has  units  of  inverse  length,  the  same  as  the  wave  number 
The  units  are  usually  called  nepers  per  unit  length. 

There  are  a  number  of  conventions  for  what  is  meant 
by  a  unit  source.  The  one  used  here  will  define  a  source  of 
strength  Q  as  one  in  which  the  pressure  p(_x)  near  the  source 
at  >t q  is  given  by 

p(x)  =  Q  |x  -  xj  ^  e^O  -^0  ^  (13) 

where  kp=  u>/c(x^).  The  pressure  due  to  the  discrete  modes 
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is  then  given  by: 


"  (U 

p(x)  =  17tQp(z)  £  u^( Zq)u^( z)Hq  '(^r) 
n=l 


(14) 


and  the  propagation  loss  relative  to  unit  distance, 
defined  by 

PL  =  -10  log10  I  p<x_)  / p(211)  1  2  (  1  5  ) 

where  U  -j  -_Xq  |  =  1,  is  given  by 

N  rn 

PL  =  -10  log1Q  | Trp ( z)  Z  un(zQ)un(z)H0  '(k^r)!2  (16) 

n=l 


This  formula  is  valid  for  all  depths  z  [8],  although  most 

authors  say  it  is  valid  only  in  layer  one. 

For  knr  large  (and  no  attenuation  in  the  bottom), 

the  expression  for  the  intensity  Ipl^  can  be  written  as, 

N 

|  p(x)  |  2  =  |  C  [  2  E  |  u  (  zn)  u  (  z)  |  2 (2/irk  r) 
ii  1  n  0  n  '  " 


n=l 


+  |  C  | 2  E  u  (z.)u  (z)u  (z  )u  (z)  - Yn 

n*n  n  0  n  n>  0  m  k  )l/2 


i(k  -k  )r 

e  n  ra  (17) 


n  m 


where  C  =  iirOp(z),  and  the  Hankel  function  has  been  replaced 
by  its  asymptotic  form 


i(k  r-rr/4) 

Hn(k  r)  ~  ( 2/irk  r)  '  e 
On  n 


(18) 


Notice  that  only  the  second  term  of  (17)  contains  the 
oscillating  trigonometric  terms.  There  are  several  reasons 
for  expecting  this  term  to  be  small.  First,  if  the  r  term 
in  the  denominator  is  removed,  the  average  value  of  the  term 
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(as  a  function  of  range)  is  zero.  Moreover,  as  a  function 
of  z,  the  mode  functions  should  appear  more  or  less  with 
random  amplitudes  and  signs.  Thus  the  main  features  of 
propagation  loss  are  contained  in  the  first  term  only.  In 
addition,  any  source  will  have  some  frequency  spread,  which 
will  smear  out  the  fluctuations  in  the  second  term.  If  the 
intensity  is  calculated  using  only  the  first  term,  a  very 
smooth  function  is  obtained,  from  which  a  more  meaningful 
comparison  of  various  source  and  receiver  depths  and  various 
ranges  can  be  made. 

This  intensity  gives  rise  to  what  is  called  the 
incoherent  propagation  loss: 

n  m 

IPL  =  -10  log  n{  (ttp(z))2  T.  ju  (z  )u  (z)H:.  ;(k  r)|2}  (1°) 

iu  .nun  u  n 

n=l 

The  incoherent  propagation  loss  should  be  used  cautiously 
for  deep  water  since  the  regions  of  high  and  low  intensity 
are  averaged  out.  In  shallow  water  it  is  a  useful  quantity 
and  for  broadband  sources  more  meaningful  than  the  coherent 
propagation  loss  PL. 

Since  the  acoustic  energy  is  trapped  by  the 
waveguide,  the  propagation  loss  as  a  function  of  range  tends 
to  behave  as  cylindrical  spreading.  It  is  sometimes  convenient 
to  remove  this  geometrical  term  by  subtracting  10  log  (r) 
from  the  propagation  loss.  Another  geometrical  factor  which 
can  be  subtracted  is  the  water  depth  correction,  10  log(h). 

This  corresponds  to  the  spreading  of  the  acoustic  energy 
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throughout  the  water  column,  and  is  closely  related  to  the 
spherical  spreading  of  the  acoustic  energy  near  the  point 
source.  The  remainder  term  contains  the  losses  due  to  the 
depth  dependence  of  the  propagation  loss  and  losses  due  to 
absorption  by  the  bottom.  If  there  is  no  attenuation  in  the 
bottom  layer,  the  remainder  term  will  be  constant  with 
range.  In  the  printout  for  the  program  (Appendix  C)  the 
incoherent  propagation  loss  with  cylindrical  spreading  and 
water  depth  correction  removed  appears  under  the  heading 
INCOH-GEOM. 

The  excitation  of  a  given  mode  n  by  a  source  at 
depth  Zq  will  be  defined  to  be 

En  =  Un(20)  (20) 

This  is  the  coefficient  of  the  n-th  term  in  the  sum  for  the 
pressure,  see  Eq . ( 1 4  )  . 


2.2  Method  of  Solution 

To  obtain  the  eigenvalues  kn  of  the  mode  equation 
(1),  the  transcendental  equation  (7)  must  be  solved.  To 
facilitate  the  solution  two  dimensionless  quantities  are 
def ined : 

y  =  h[w2/c2  -  k2]^2  (21) 

7n  1  n 

and 

b  =  (hw/c^)  [1  -  c^/c2]1^2  (22) 
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(23) 


(214) 


A  Newton-Raphson  iteration  method  is  used  to  solve 
Eq  .  (214)  for  yp.  To  first  get  a  feeling  for  the  results,  a 
graphical  solution  to  Eq.(23)  can  be  obtained.  This  is 
illustrated  in  Figure  2  for  the  case  of  three  modes.  Notice 
that  the  roots  yn  must  occur  on  the  negative  branch  of  the 
tangent  function  or 

(n-4)lr<y<mr  (25) 

1  n 
for  n= 1 , 2, . . . , N,  where 

N  -  [[  \  +  £  ]]  (26) 

and  the  notation  ffx]J  means  the  greatest  integer  less  than 
or  equal  to  x.  In  terms  of  the  environmental  input 
parameters  Eq.(26)  becomes 

N  =  [[  |  +  (2hf/cx)  (1  -  c \lc\)112  ]]  (27) 

Also,  since  the  right  hand  side  of  Eq.(22)  is  monotonicall y 
decreasing,  it  follows  that  successive  values  of  yn  differ 
by  less  than  n.  So  yn  is  constrained  by 

(n  -  4)  ti  <  y  <  min  (nn,  b,  y  ,  +  n)  (28) 

2  n  n-1 

For  the  Newton -Raphson  method  we  use  the  function 


F(y)  defined  by 


F(y)  =  y  +  tan  L  [  — —— — — -y-  ]  (29) 

pi  (b2-y2) 1/2 

The  derivative  of  F  is  then  given  by: 

( p  /p  ) 

F’(y)  =  1  + - rj- - — - - (30) 

(b2-y2)  '  {l+[ (p2/p1)2-1] (y/b)2} 

and  an  improved  estimate  of  yn  is  obtained  from  the 

Newton-Raphson  iteration  formula: 


(1+1) 

y„ 


F(y^i))/F,(y^i)) 


(31  ) 


In  the  numerical  procedure  we  must  ensure  that  at 

each  iteration  the  new  value  for  yn  remains  in  the  range 


Once 


(n  -  -|)  ti  <  <  rain  {mi,  b)  (32) 

has  been  determined,  the  wave  number  is  given  by 


k  =  [u2/c2  -  y2/h2]^^  (33) 

n  I  n 

The  mode  functions  and  propagation  loss,  etc.,  can  then  be 
calculated  from  the  equations  of  Section  2.1. 
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3.  PROGRAM  DESCRIPTION 
3.1  Program  Structure 

The  program  N0RM2L  is  written  in  modular  form  for 
ease  of  expansion  and  modification.  An  attempt  has  been  made 
to  remain  close  to  standard  Fortran,  although  for  input  and 
output  some  straying  did  occur.  In  particular,  for  input  of 
data,  free  field  formatting  was  used;  and  for  printing  of 
text,  quotes  were  often  used  instead  of  the  Hollerith 
specification.  A  simplified  flowchart  of  N0RM2L  is  given  in 
Figure  3,  and  the  interaction  of  the  subroutines  is  depicted 
in  the  tree  structure  of  Figure  4.  Appendix  A  is  a  User's 
Guide  for  program  N0RM2L  with  details  on  the  hardware- 
software  environment,  input/output,  operational  instructions, 
restrictions  and  limitations,  functions  and  subroutines, 
accuracy,  memory  requirements,  and  execution  times. 

Below  is  a  list  of  the  subroutines  and  COMMON  blocks 
arid  a  one  line  description  of  their  functions.  A  short  but 
more  detailed  description  of  t.he  functions  and  subroutines 
appears  in  Appendix  A. 6.  The  COMMON  blocks  are  discussed 
further  in  Section  3*F. 
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List  of  Modules  and  COMMON  Blocks 


(a)  Modules  included  in  N0RM2L 


N0RM2L 

NMODES 

MODE 

ATTENU 

MODFUN 

WFUNCT 

UN 

PLOSS 
EXC ITE 

PROPLO 
VPH I 

HANK01 


Main  program.  Driver  program  for  I/O 
Calculates  modes  by  calling  MODE 
Calculates  eigenvalue  for  a  single  mode. 
Calculates  attenuation  for  the  modes. 

Prints  modal  functions  at  12  depths 
Calculates  mode  functions  at  specified  depths 
Evaluates  a  mode  function  at  a  specified 
depth 

Propagation  loss  I/O  routine 

Calculates  mode  excitations  for  a  specified 

source  depth 

Calculates  propagation  loss  in  dB//unit  length 
Calculates  velocity  potential  at  the  specified 
receiver  depth 

Evaluates  the  asymptotic  Hankel  function 


(b)  DEC-20  subroutines  (not  included) 

DATE  Returns  current  date 

TIME  Returns  current  time 


(c)  COMMON  blocks 


/ CNTLIO/ 
/ENVIRN/ 
/EIGENF / 
/CNTLNM/ 


I/O  control.  Contains  logical  unit  numbers. 
Environmental  and  derived  quantities 
Eigenvalues  and  normalization  constants 
Control  parameters  for  calculating  modes 


COMMON  block  usage: 


COMMON  Block  Modules  that  use  the  Block 


/CNTLIO/ 

/ENVIRN/ 

/EIGENF/ 

/CNTLNM/ 


N0RM2L, NMODES, MODES, PLOSS 
N0RM2L , ATTENU , MODF UN , WF UNCT , PLOSS , VPH I 
N0RM2L, ATTENU, WFUNCT, PLOSS, VPH I 
N0RM2L, NMODES, MODES 


i 

f 


3.2  Important  variables 

3.2.1  Variables  in  COMMON  blocks 


COMMON 

/ENVIRON/ 

Quantity 

Variable 

name 

Default  Description 

value 

h 

H 

100. 

thickness  of  first  layer 

C1 

C  1 

1500. 

sound  speed  in  first  layer 

Pi 

RHO  1 

1 . 0 

density  of  first  layer 

c2 

C2 

1800. 

sound  speed  in  bottom  layer 

P2/  P-| 

RHORAT 

2.0 

ratio  of  densities 

e 

ATTN 

CM 

O 

II 

CD 

attenuation  in  bottom  layer 
at  given  frequency 

oj  =  2  it  f 

OMEGA 

f=100 

.  angular  frequency 

k  1  =  00/ c  1 

AK1 

wave  number  in  first  layer 

k  2=  C0/C2 

AK2 

wave  number  in  bottom  layer 

COMMON 

/CNTLNM/ 

Variable 

name 

Default 

value 

Description 

MAX 

20 

maximum 
in  the 

number  of  iterations  allowed 
Newton-Raphson  procedure 

TOL 

1 . 0E-8 

the  tolerance  desired  in  computing  yn 

*NMIN 

1 

first  mode  to  be  calculated 

NMA  X 

50 

maximum  number  of  modes  that  can  be 
calculated 

*VM  IN 

C1 

minimum 

phase  velocity  of  modes 

*VMA  X 

c2 

max imum 

phase  velocity  of  modes 

*  Treated 
expansion 

as  constants 

II 

in  this 

version.  "Reserved  for  future 
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COMMON  /EIGENF / 


Quantity 

Variable 

name 

Description 

N 

NM 

The  number  of  modes  actually  calculated 

kn 

AKM ( 50  ) 

Array  containing  wave  numbers  of  the  modes 

An 

ANM ( 50 ) 

Array  containing  normalizations  of  the 
mode  functions 

«n 

ATN (50) 

Array  containing  imaginary  part  of  the 
wave  numbers 

COMMON 

/CNTLIO/ 

Logical 

Current  implementation 

File  type 

var iable 

Logical  file 

Device 

Input 

II 

5 

Terminal 

Output 

IT 

5 

Terminal 

Output 

10 

3 

Line  pri 

There  are  three  logical  files  for  input/output 
although  in  the  current  version  two  of  these  are  identical. 
Variable  names  are  used  rather  than  numbers  to  allow  for 
compatibility  with  other  computers  and  other  input/output 
devices.  For  example,  on  DEC'S  PDP— 1 1 3^  computer  the  logical 
file  numbers  are  7,  7,  and  6  respect i vel  y .  Only  one  data 
statement  needs  to  be  changed.  It  may  even  be  possible  to 
run  the  program  in  batch  mode  by  letting  II  correspond  to 
the  card  reader,  IT  to  some  arbitrary  file,  and  10  to  the 
line  pr inter  . 
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3-2.2  Other  important  variables 


Quantity 

Variable 

name 

Default 
val  ue 

Description 

c  -  /  c  i 

CRAT 

1 . 2 

ratio  of  sound  speeds 

2 

RH02 

P2/P i =2 • 

density  in  bottom  layer 

a 

ALPRAT 

0.2 

attenuation  coefficient  in 
bottom  layer  (dB/length 
//  1kHz) 

f 

FREQ 

100. 

frequency 

ZC 

ZS 

h/2 

source  depth 

z 

ZR 

h/2 

receiver  depth 

r 

R 

range  of  receiver 

RMIN 

10000. 

minimum  range  for  propagation 
loss  (PL)  calculations 

DELR 

10000. 

range  increment  for  PL 
calculations 

NRNG 

10 

number  of  range  increments 
for  PL  calculations 

vn 

VPHASE 

phase  velocity  for  mode  n 

PL 

PL 

propagation  loss  in  dB 
(relative  to  unit  distance) 

I  PL 

I  PL 

incoherent  PL  in  dB 

r  .  I PL 

RCOHI 

IPL  with  cylindrical 
spreading  removed 

n 

pr 

Greek  pi  (3.  1«M59265  ) 

h  r  i 

HOMO 

useful  dimensionless 
quantity 

b  ,  b‘ 

B,  BSQ 

useful  quantity  for 
determining  eigenvalues 
(see  Eq . 22  ) 

17 


3-3  Testing  and  Verification 


To  verify  that  the  computer  program  was  giving 
accurate  answers,  its  output  was  compared  with  two  other 
sources:  u  HP-67  programmable  hand  calculator  program  and  a 
normal  mode  computer  program  by  Bartberger  and  Ackler  [9] 

(to  be  called  BANORM  in  this  report).  The  HP-67  program 
computes  only  the  eigenvalues  using  an  algorithm  similar  to 
the  one  described  in  Section  2.2 ,  but  the  calculator  carries 
about  10  significant  digits,  compared  with  about  8  for  the 
I'Kf-.’U  computer  ;  thus  its  solution  will  be  regarded  as 
exact.  The  program  BANORM  is  a  general  normal  mode  program 
which  calculates  eigenvalues,  mode  functions  and  propagation 
loss  for  more  complicated  sound  speed  and  density  profiles 
than  the  Pekeris  model.  It  was  chosen  for  comparison 
because  it  was  available,  and  because  the  C DC  6400  computer 
on  which  it  was  implemented  had  a  precision  of  14 
significant  digits.  The  program  used  a  numerical 
integration  of  the  differential  equation,  so  its  accuracy 
was  not.  known  beforehand.  However,  as  will  be  seen,  its 
eigenvalues  seemed  correct  so  its  mode  functions  and 
propagation  loss  were  assumed  to  also  be  correct.  In  the 
discussion  below,  the  eigenvalues  from  the  HP-67  program  and 
the  eigenvalues,  mode  functions  and  propagation  loss  from 
BANORM  are  compared  with  the  cor  respond ing  values  from 
NORM,  1. . 

Several  modifications  to  the  normal  mode  program 
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BANORM  were  necessary  before  the  results  could  be  compared. 
First,  the  correction  to  the  sound  speed  profile  due  to  the 
earth's  curvature  had  to  be  removed.  Removal  of  the 
curvature  caused  the  WKB  approximation  procedure  to  fail,  so 
an  additional  modification  was  made.  Another  correction 
removed  was  the  additional  propagation  loss  due  to  volume 
absorption  in  the  water  column.  One  difference  which  could 
not  be  removed  was  the  treatment  of  absorption  in  the  bottom 
layer.  Program  BANORM  treats  the  absorption  "exactly"  by 
using  an  imaginary  value  for  the  wave  number,  whereas  in 
N0RM2L  the  absorption  is  treated  as  a  perturbation.  The  two 
programs  should  give  identical  results  for  the  case  where 
there  is  no  attenuation  in  the  bottom.  To  facilitate 
comparison  of  results  BANORM  was  modified  to  accept  input 
and  print  output  in  metric  units,  although  the  internal 
calculations  were  performed,  as  before,  in  yards.  The  values 
ir.  the  tables  have  been  modified  using  the  exact  conversion 
factors  as 

k  >  k  / (0.9 144) 
n  n 

u  (z)  (-l)n+1(0.9144)“i/2u  (z) 

n  n 

PL  >  PL  +  logi0(0.9144) 

Table  I  shows  the  results  from  the  HP-67  two  layer 
program,  with  the  default  values  for  the  model.  These  values 
should  have  about  9  significant  figures  accuracy  and  were 
user,  to  check  output  from  the  computer  program.  A  comparison 


with  the  results  shown  in  Appendix  B  shows  excellent 
agreement,  with  eight  digits  accuracy  in  kn. 


Table  II  compares  the  eigenvalues  kn  calculated  by 
NCRM2L  and  BANORM  for  the  cases  of  no  attenuation  and 
attenuation  equal  to  P.2  dR/m//1  kHz.  The  results  of  the 
HP-67  calculation  are  repeated  for  easy  comparison.  For  no 
attenuation  all  three  methods  agree  to  about  eight 
significant  digits,  and  the  HP  and  CDC  results  agree  to  the 
full  nine.  Since  the  HP  program  did  not  calculate  the  mode 
functions  or  propagation  loss,  BANORM  was  used  as  the 
reference  for  judging  the  accuracy  of  N0RM2L.  When  the 
bottom  absorption  is  non-zero,  the  exact  solution  shows  a 
shift  in  the  real  part  of  the  eigenvalue  as  well  as  the 
introduction  of  an  imaginary  part.  The  shift  in  the  real 
part  is  less  than  the  imaginary  part,  however.  In  the 
example  the  real  parts  of  kn  agree  to  about  six  figures 
accuracy  for  the  exact  and  perturbation  calculations.  The 
imaginary  parts  agree  to  about  three  figures  accuracy,  or 
about  the  seventh  decimal  place.  Notice  that  the  imaginary 
part,  though  small,  increases  with  increasing  mode  number  - 
a  factor  of  about  60  in  the  seven  modes. 

The  mode  functions  are  compared  at  two  depths  in 
Table  III.  The  answers  agree  to  about  six  decimal  places. 
Presumably,  i- A  NORM  is  the  more  accurate  because  of  the 
erect  er  accuracy  1  r.  the  wave  numbers  kn.  As  is  quite  common 
in  numerical  approximations,  the  eigenfunctions  are  less 


accurate  than  the  eigenvalues.  When  attenuation  is 
considered,  the  mode  functions  agree  to  only  four  digits. 

Table  IV  compares  the  propagation  loss  calculated  by 
N0RM2L  and  BANORM  for  source  and  receiver  at  mid-depth,  with 
bottom  absorption  both  present  and  absent.  The  values  for 
the  average  or  incoherent  propagation  loss  (I PL)  are 
considerably  more  in  agreement  than  those  for  the  coherent 
propagation  loss  (PL)  -  five  digits  compared  with  four  for 
no  ebsorption,  and  four  digits  compared  with  three  digits 
for  attenuation  coefficient  a=0.2.  For  the  IPL  there  seems 
to  be  a  consistent  bias  with  BANORM  giving  values  low  by 
0.0002  dB  for  a=0  and  by  C.002  for  a=0.2.  At  present  no 
explanation  can  be  given  for  this  bias. 

In  summary,  for  no  bottom  absorption  the  wave 
numbers  as  calculated  by  N0RM2L  are  probably  accurate  to 
eight  significant  digits,  the  mode  functions  accurate  to 
about  six  digits,  and  the  propagation  loss  to  five  digits. 
When  absorption  is  included  as  a  perturbation,  the  accuracy 
of  the  results  is  degraded,  but  the  propagation  loss  remains 
accurate  to  better  than  0. 1  dB. 
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4.  CONCLUSIONS  AND  RECOMMENDATIONS 


NORM2L  was  written  to  provide  an  interactive 
computer  program  for  calculating  the  discrete  normal  modes 
and  acoustic  propagation  loss  for  the  Pekeris  model  of  the 
ocean.  Although  the  model  is  simple,  it  contains  all  the 
essential  features  of  normal  mode  theory,  and  for  many 
shallow  water  environments  is  a  reasonable  approximation  to 
the  actual  physical  situation.  The  simple  model  has  the 
advantage  that  the  extensive  calculations  of  a  general 
normal  mode  program  are  not  necessary  and  an  interactive 
program  can  be  used.  The  interactive  nature  of  the  program 
makes  it  useful  for  sensitivity  studies  since  the  user  can 
see  the  results  immediately,  change  the  input  parameters  and 
look  at  a  new  set  of  results.  The  user  can  then  quickly  get 
a  feeling  for  the  importance  of  the  physical  parameters. 

In  practice  the  program  has  proved  to  be  very  useful 
as  a  learning  tool  for  normal  mode  calculations  and  for 
sensitivity  studies.  In  view  of  the  utility  of  the 
program  it  is  recommended  that  certain  extensions  and 
improvements  be  made  to  the  program.  Certain  improvements 
could  be  made  to  the  present  program  to  provide  more  options 
and  improve  output  to  include  graphing  of  mode  functions  and 
propagation  loss.  These  are  listed  in  Appendix  E.  Some 
extensions  could  be  made  to  the  model  to  improve  its 
usefulness  for  sensitivity  studies  by  allowing  for  a  more 
realistic  environment.  For  example,  propagation  loss  will 
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be  affected  by  the  presence  of  shear  waves  in  the  bottom  or 
scattering  from  a  rough  surface  or  bottom.  These  could  be 
included  as  perturbation  effects  to  be  added  to  the  present 
program,  or  treated  more  rigorously  in  an  extended  program. 
If  the  acoustic  field  is  needed  near  the  source,  the 
continuous  modes  are  important  and  should  be  calculated. 
Range  dependence  could  be  included  in  the  model  to  allow  for 
a  sloping  bottom  or  changing  bottom  properties.  This  could 
be  incorporated  in  an  approximate  manner  in  the  adiabatic 
approximation  (no  mode  coupling),  or  more  exactly  by 
including  mode  coupling  effects. 
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Table  I.  Eigenvalues  obtained  using  the  Newton-Raphson 
iteration  procedure 


mode 

iterations  yn 

kn 

vn 

1 

10 

2 . 894719966 

0.417877606 

1503.594646 

2 

9 

5.805255120 

0.414836757 

1514.616340 

3 

9 

8 . 740763033 

0. 409657833 

1533.764231 

4 

8 

1  1.70219503 

0.402200817 

1562.201031 

5 

7 

14. 68458504 

0.392295718 

1601 . 645144 

6 

7 

17. 67884983 

0. 37974394 1 

1654. 584742 

7 

8 

20.66772252 

0.364340718 

1724.535581 

The 

results 

shown  are  for  the 

HP-67  program 

with  the  default 

parameters  h=100,  c-|  =  15C0,  p i  =  1  ,  C2/c-|  =  1.2,  P2/  P-|=2.0,  and 

f=100.  The  results  were  obtained  by  a  Newton-Raphson 
procedure  which  iterated  until  successive  approximations  to 
yn  differed  by  less  than  1.0E-9.  The  initial  approximation 
was  yn  =  (n-1/2)ir  . 


Table  IT.  Comparison  of  Wave  Numbers 


Mode 

Exact 

N0RM2L 

Bartberger-Ackler 

kn 

kn 

6 

n 

kn(a=0) 

kn(a=. 2) 

<$  n(a=.  2) 

1 

.417877606 

.41787760 

.23081 (-5) 

.417877606 

.417877551 

.23066 (-5) 

2 

.414836757 

.41483675 

.85045 (-4) 

.414836757 

.414836565 

.84998 (-4) 

3 

.409657833 

.40965783 

.17149 (-4) 

.409657833 

.409657434 

.17141  (-4) 

4 

.402200817 

.40220081 

.27377 (-4) 

.402200817 

.402200241 

.27364 (-4) 

5 

.392295718 

. 39229571 

.39711 (-4) 

.392295718 

.392294853 

.39688 (-4) 

6 

.379743941 

.37974394 

.  57252 (-4) 

.379743 941 

.379742482 

.57202 (-4) 

7 

.364340718 

.36434072 

.94016 (-4) 

.364340718 

.364336855 

.93761 (-4) 

Ttie  default  parameters  given  with  Table  I  have  been  used  for  the 
calculations.  The  exact  values  are  from  the  HP-67  program. 
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Table  III 

.  Comparison 

of  Mode  Functions 

At  depth 

50 

mode 

N0RM2L 

BAN0RM ( a  =  0  ) 

BANORM ( a=0 

1 

. 1348502 

.  1348504 

. 1348530 

2 

. 32284 1 9 ( - 1  ) 

. 3228438 ( - 1  ) 

.  3227565  ( 

3 

. 1290802 

1290816 

1290872 

4 

.  576 1266 (-1  ) 

- -  576  1 346  (- 1  ) 

5760 1 86 ( 

5 

. 1203002 

. 1203022 

. 1203128 

6 

. 7630017 (-1  ) 

.  76301  39  (-1  ) 

. 7628658 ( 

7 

. 1085481 

1085494 

1085981 

At  depth 

25 

mode 

N0RM2L 

BANORM ( a=0  ) 

BANORM ( a= 

1 

. 89975 1 5 ( - 1 ) 

.  8997537 ( - 1  ) 

.8997854 ( 

2 

. 1354220 

. 1354228 

.  1  354246 

3 

. 1 119586 

. 1 119598 

.1119571 

4 

.  2949  t  92 (- 7  ) 

. 2949233 (- 1 ) 

. 2948582 

5 

.  6969622 ( - 1  ) 

6969738 (- 1  ) 

-.6970512 

b 

.  1  322287 

1322308 

1322377 

7 

. 1236168 

1236184 

-. 1236284 

.2) 

-1  ) 

-1  ) 
-1  ) 


0.2) 


-1  ) 


The  default  parameters  given  with  Table  I  were  used  for 
calculations . 


Table  IV. 


Comparison  of  Propagation  Loss. 


a  =  0 


NORM  Li.  RANORM 


r 

PL 

IPL 

PL 

I  PL 

10000 

57. 4530 

58.0333 

57.4524 

58.0331 

20000 

62. 1264 

61.04  36 

62. 1259 

61 . 0434 

30000 

62. 9224 

62.8047 

62. 9221 

62.  8043 

40000 

64 . 3649 

64.0539 

64. 3650 

64 . 0537 

50000 

81.6223 

65.0230 

81.6213 

65.0228 

60000 

66. 5136 

65.8148 

66.5127 

65 • 8146 

70000 

73- 5068 

66. 4842 

7^.5066 

66.4841 

80000 

64.2857 

67. 0642 

64.2868 

67. 0640 

90000 

72.6916 

67.5757 

72.6934 

67.5755 

1000C0 

64 . 7800 

68.0333 

64 . 7807 

68.0331 

a  =  0 . 


NORM  2L  BANORM 


r 

PL 

IPL 

PL 

IPL 

10000 

58.3779 

60.0250 

58. 4343 

60.0225 

20000 

64. 9796 

64. 2338 

65. 1 382 

64.2318 

30000 

65. 4292 

66.8262 

65.4717 

66. 8244 

40000 

70.6964 

68.7062 

70.6569 

68.7044 

50000 

76. 3537 

70. 1800 

76.4301 

70. 1782 

60000 

75. 9126 

71 . 3920 

76. 0060 

71 . 3903 

70000 

74.7604 

72.4229 

74.7096 

72.421  1 

80000 

74.8863 

70. 3220 

74.7803 

73. 320? 

90000 

74.897  1 

74. 1222 

74.8337 

74. 1203 

100000 

74.4043 

74.8458 

74 . 5660 

74.8439 

The  standard  environmental  model  is  used;  source  and 
receiver  are  at  mid-depth.  Pi.  is  the  coherent  propagation 
loss,  IPL.  is  the  incoherent  propagation  loss,  and  a  is  the 
attenuation  coefficient  . 
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DEPTH 


RANGE  r 

PRESSURE  RELEASE  SURFACE  ^ 

z  --  zs  X  SOURCE 

z  --  zr  RECEIVER* 

Pt  ,<=, 

(a)  The  layer  structure:  densities  p-j  and  pp,  sound 
speeds  c-j  and  €£,  and  bottom  absorption  ap  are  all  constants, 
and  c^  <  C2<  The  source  and  receiver  are  at  depths  zs  and  zr 
respectively  and  separated  by  horizontal  distance  r. 


RANGE  0 

z  -  0  — - •— 


z  =  0 


SOUND  SPEED 


x 

t— 

Q. 

UJ 

Q 


Z  =  h 


(b)  The  sound  speed  profile  for  the  Pekeris  model. 


Figure  1.  The  Pekeris  model 
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Figure  2.  Graphical  Solution  for  the  Mode  Eigenvalues 

Graphical  solution  of  Eq.(23),  for  P2/p1=2.0  and 
b=10.  The  solid  line  is  the  tangent  function.  The  dashed 
line  is  the  function:  -(  p-]/  Pp )  y/ [  b^-y  ^  ^  .  There  are 

three  modes  in  this  example.  The  solutions  yn  are 


indicated . 
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A.1  Identification 

Pr of r am  :  N0RM2L 

Purpose  :  To  calculate  discrete  normal  modes  and  propagation 
loss  for  a  two  layer  fluid  acoustical  model. 

Author  :  Dale  D.  Ellis 

Installation  :  Defence  Research  Establishment  Atlantic 

P.0.  Box  1012 
Dartmouth,  N.S. 

Canada,  B2Y  3Z7 

Computer:  DEC-20 

Language:  Fortran  IV;  with  some  minor  DECsystem-20 
extensions  for  input/output. 

Date  written  :  April  1978 

Last  Modification  :  June  1978 

Present  version  :  3-0 

Security  :  Unclassified 


A. 2  Hardware-Software  Environment 

Hardware  environment: 

-  Interactive  terminal 

-  Line  printer  (optional) 

Software  environment: 

-  Interactive  operating  system 

-  FORTRAN  IV  compiler  (with  a  few  DEC-20 
extensions  for  input/output  statements) 

-  Standard  Fortran  library  functions 

-  DEC-20  subroutines  for  date  and  time. 


Note:  None  of  the  DEC-20  non-standard  items  are 
essential  to  the  program;  they  are  for  convenience  only  and 
may  be  omitted  or  replaced  by  functions  compatible  with  the 
operating  system  being  used.  The  most  difficult  item  to 
change  will  be  the  free-field  format  input  and  non-ANSI 
output,  but  these  items  are  becoming  quite  common  in  FORTRAN 
comp i 1 e  r  s  . 
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A . 3  Input/Output 
A . 3 ■ 1  Input 

Most  input  is  list  directed;  that  is  no  formatting 
is  required  [10]. 

The  exception  to  this  is  the  response  to  a 
question  which  requires  a  YES  or  NO  answer 
In  this  case  no  spaces  are  allowed.  The 
program  checks  only  to  see  if  the  first 
character  is  a  Y,  otherwise  NO  is  assumed. 

If  no  change  is  desired  to  a  data  item  (other  than  Y/N)  a 

null  entry  (not  a  zero)  will  allow  its  value  to  remain 

unchanged.  On  the  DEC-20  a  comma  (,)  is  used  to  signify  a 

null  entry.  (Note:  a  carriage  return  by  itself  is  ignored.) 

The  data  items  3nd  their  default  values  are: 


Vsr i able 

Fortran  name 

Default  value 

depth  of  layer  1 

H 

100. 

sound  speed  in  layer  1 

C  1 

1500. 

density  in  layer  1 

RHO  1 

1 . 0 

ratio  of  sound  speeds 

CRAT 

1  .  2 

flayer  2  to  layer  1) 

r  . t i o  of  densities 

RHORAT 

2.  C 

attenuation  coefficient 

ALPRAT 

0.2 

for  layer  2 

frequency 

FREw 

100. 

source  depth 

•?  q 

H/2. 

receiver  depth 

ZR 

H/2 . 

minimum  range 

RMIN 

10000. 

range  increment 

DEL'ri 

10000. 

number  of  ranges 

NRNG 

10 

The  above  values  are  typical  units  for  isovelocity  shallow 
water  (layer  1)  over  a  hard  bottom  (layer  2).  The  depths 
and  ranges  are  in  meters,  the  speeds  in  meters  per  second 
and  the  frequency  in  Hz.  The  attenuation  coefficient  is  in 
>iP/m  at  a  frequency  of  1  kHz.  The  attenuation  is  assumed  to 
vary  linearly  with  frequency.  Units  need  not  be  metric  but 


should  be  consistent.  In  particular,  the  depths  and  ranges 
should  be  in  the  same  units. 

Some  checking  of  the  input  parameters  is  performed 
but  it  is  by  no  means  exhaustive.  All  input  parameters 
should  be  positive,  except  ALPRAT  which  may  be  zero.  If 
! RHC 1-1 . 0 ] >0 . 03  a  confirmation  is  required.  If  a  negative 
value  is  entered  for  the  frequency,  no  calculations  are 
performed.  This  enables  the  user  who  notices  a  previous 
incorrect  entry  to  begin  again.  Since  list  directed  input 
is  used,  only  those  values  that  need  altering  need  be 
entered.  For  the  others  a  null  entry  (i.e.  a  comma,  not  a 
zero)  and  carriage  return  will  suffice. 

YES/NO  entries  are  required  for  control  of 
printout  and  calculation  of  mode  functions,  mode 
excitations  and  propagation  loss. 

In  the  propagation  loss  subroutine  the  default 
values  are  reset  each  time  it  is  called.  The  parameters  and 
their  default  values  were  listed  earlier.  Other  values  may 
be  selected  and  the  propagation  loss  subroutine  may  be 
called  repeatedly  with  different  values  for  the  source  and 
receiver  depths  and  range  selections. 

A . 3 • 2  Output 

Output  consists  of  a  summary  of  the  input 
parameters  for  the  calculation,  the  date  and  time,  and 
results  associated  with  the  discrete  normal  modes.  In  the 
summary  MODE  is  the  mode  number,  KN  is  the  mode  wave  number, 
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NORM  is  the  normalization  of  the  mode  functions,  IM(KN)  is 
the  imaginary  part,  of  the  wave  number  and  VPHASE  is  the 
phase  speed.  The  output  is  normally  displayed  on  the 
terminal.  However,  if  YES  is  specified  in  answer  to: 

'DO  YOU  WANT  A  COPY  OF  THIS  SUMMARY  ON  THE  LINE  PRINTER  ?' 
then  all  further  output  from  the  calculations  of  the  mode 
functions,  propagation  loss  and  mode  excitations  is  printed 
on  the  line  printer.  Only  the  prompting  questions  continue 
to  be  printed  on  the  terminal. 

The  mode  functions,  if  desired,  are  calculated  at  12 
depths  equally  spaced  by  H/10,  where  H  is  the  thickness  of 
layer  one.  The  first  line  of  output  lists  the  depths.  The 
subsequent  lines  list  the  mode  functions  at  the 
corresponding  depths.  Note  that,  due  to  the  change  in 
density,  the  mode  functions  are  discontinuous  at  the 
interface  between  the  two  layers. 

The  mode  excitations  are  simply  the  normalized  mode 
functions  evaluated  at  the  source  depth.  They  may 
optionally  be  printed  with  the  propagation  loss  calculation. 

At  each  range,  coherent  propagation  loss,  incoherent 
propagation  loss,  and  incoherent  propagation  loss  with 
cylindrical  and  (approximate)  near-field  spherical  spreading 
removed  are  calculated  and  printed  in  tabular  form  under  the 
headings  RANGE,  COH,  INCOH,  and  INCOH-GEOM.  The  propagation 
loss  subroutine  can  be  called  repeatedly  with  different 
values  for  the  source  depth,  receiver  depth,  and  range 
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points.  Note  that  each  time  the  subroutine  is  called  the 
values  revert  to  the  default  values  listed  in  the  previous 
section . 

Sample  printout  from  a  terminal  is  shown  in  Appendix 
B.  The  underlined  characters  are  entered  by  the  user.  Data 
entries  are  terminated  by  pressing  the  "return"  key;  this  is 
indicated  by  <ret>  on  the  output. 

Appendix  C  shows  the  line  printer  output 
corresponding  to  the  input  shown  in  Appendix  B.  If  the  line 
printer  option  is  not  specified,  this  printout  will  appear 
at  the  terminal. 

A. 4  Operational  Instructions 

Initiate: 

This  is  an  interactive  program;  no  special 
instructions  are  necessary.  For  example  on  the  DEC-20,  if 
the  program  has  been  SAVEd  as  file  N0RM2L.EXE  on  the 
directory  <ELLIS.NM>,  then  the  command 

@RUN  <ELLIS. NM>N0RM2L 

will  cause  execution  to  begin.  No  system  messages  appear. 

In  the  example  of  Appendix  B  the  compiled  program  exists  on 
file  N0RM2L  of  the  user's  directory,  and  the  command 

^EXECUTE  N0RM2L 


causes  the  relocatable  files  to  be  linked  with  the  system 
subroutines  and  loaded  for  execution.  There  are  two  lines 
of  system  messages.  If  only  the  source  file  N0RM2L.F0R  had 


been  on  disk,  then  the  FORTRAN  compiler  would  have  been 
called  automatically  and  each  subroutine  name  would  have 
been  typed  as  ir  was  compiled. 

Execute : 

Enter  data  in  response  to  questions  typed  on  the 
terminal.  See  Appenxix  A. 3  for  more  details.  As  many  data 
sets  as  desired  may  be  executed. 

Input/output  data: 

Both  input  and  output  appear  on  the  terminal. 
Optionally,  some  output  may  be  directed  to  the  line 
printer.  See  Appendices  B  and  C  for  an  example;  responses 
typed  by  the  user  are  underlined. 

Error  procedures: 

Error  checking  is  not  exhaustive  although  those 
errors  that  are  trapped  allow  the  user  to  change  user  input. 
Some  warning  messages  are  issued  -  see  below  under 
"Messages".  The  DEC-20  operating  system  traps  overflow, 
square  roots  of  negative  numbers,  etc.,  and  issues  warning 
messages  for  the  first  two  errors  detected.  Any  system 
error  messages  indicate  that  the  results  may  be  unreliable. 

Interrupt : 

No  interrupt  procedures  have  been  coded  into  the 
program.  However,  on  the  DEC-20  typing  ’control-C'  once  or 
twice  will  cause  program  execution  to  stop.  The  program  can 
be  re-run  with  the  RUN  or  EXECUTE  command  as  described  in 
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the  "Initiate"  paragraph  above.  The  START  or  CONTINUE 
commands  can  be  used  to  restart  the  program  without 
reloading,  but  these  should  not  be  used  if  an  error  has 
occurred . 


Terminate : 

Program  execution  can  be  terminated  by  typing  NO 
in  response  to  the  question 
’DO  YOU  WISH  TO  RUN  ANOTHER  DATA  SET  ?  (Y/N)’. 


Messages : 

The  following  warning  messages  are  issued: 

Message  Explanation  Response 

ARE  YOU  SURE  IRH01-1.I  >  0.03  Answer  YES  if  RH01 

ABOUT  RH01?  (Y/N)  Usually  the  density  was  entered  correctly; 

is  entered  as  the  a  NO  will  cause  the 
specific  gravity  previous  question 

1.0  or  1.026.  This  to  be  repeated, 
message  is  to  prevent 
a  zero  or  ridiculous 
value  being  entered 
by  mistake. 

NEGATIVE  FREQUENCY.  A  negative  number,  User  is  asked  if  he 
EXECUTION  either  deliberately  wants  to  run  another 

SUPPRESSED.  or  by  mistake,  was  data  set. 

entered  for  the 
frequency.  No  modes 
are  calculated. 

***  WARNING  ***  In  the  present  None.  User  must 

THERE  ARE  [ NM ]  version  the  array  decide  whether  the 

MODES.  ONLY  [NMAX]  sizes  have  been  neglected  modes  will 

MODES  HAVE  BEEN  dimensioned  for  50  affect  his  results. 
CALCULATED.  modes.  Program 

calculates  the  NMAX 
modes  and  continues. 


NO  TRAPPED  MODES.  There  are  no  User  is  asked  about 

discrete  modes  for  the  usual  line 
the  current  set  of  printer  summary  and 
environmental  whether  to  run 

parameters.  Further  another  data  set. 
calculations  are 
skipped . 

The  Newton-Raphson  Indicates  TOL  has 

procedure  for  mode  been  chosen  too 

M  has  failed  to  small,  MAX  not  large 

converge  to  the  enough,  or  a  bug  in 

desired  tolerance  the  program.  See 
after  MAX  iterations  section  3-3  for 
The  difference  details  on  how  to 

between  the  last  modify  TOL  or  MAX. 

two  values  of  YN  is 
DIF.  Answers  may  be 
inaccurate . 


A. 5  Restrictions  and  Limitations 

The  limitations  of  the  model  are  discussed  in 
Section  2.1.  This  section  refers  to  the  restrictions  and 
limitations  of  the  program. 

All  input  parameters  must  be  positive;  except 
ALPRAT  which  may  equal  zero. 

The  sound  speed  in  the  bottom  must  be  greater  than 
the  sound  speed  in  the  water  column;  otherwise,  no  modes 
will  be  trapped. 

The  present  version  of  the  program  allows  for  a 
maximum  of  50  modes.  This  can  be  modified  by  changing  a  few 
DIMENSION  statements. 

A  small  value  for  the  range  variable  may  cause  the 
Hankel  function  to  be  inaccurate.  No  warning  is  given.  The 
error  would  occur  at  short  ranges  where  the  model  is 


CONVERGENCE  NOT 
ACHIEVED  FOR  MODE 
[M]  AFTER  [MAX] 
ITERATIONS,  DIF= 
[DIF]. 
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inaccurate  anyway. 


The  program  is  written  in  single  precision.  See 
Appendix  A. 7  for  more  comments  on  accuracy. 

Printout  of  the  mode  functions  and  propagation  loss 
appears  on  the  line  printer  or  terminal,  but  not  both. 


A. 6  Functions  and  Subroutines 

A  chart  of  the  subroutine  interactions  is  given  in 
Section  3-1  •  The  program  was  designed  to  be  modular  for 
flexibility  in  expansion  and  modification;  thus  some 
duplication  and  inefficiency  results. 

The  program  has  the  following  functions  and 
subroutines.  For  any  undefined  variables  see  Section  3.2. 

1.  subroutine  NMODES  (H , HOMC , CRAT , RHO  1 , RHORAT , NM , AKM , ANM  ) 

-  given  the  input  environmental  parameters  H,  HOMC, 
CRAT,  RH01,  and  RHORAT,  calculates  the  number  of  trapped  modes 
(NM),  the  wave  numbers  of  the  modes  (array  AKM),  and  the 
normalization  factors  for  the  wave  functions  of  each  mode. 


2.  subroutine  MODE  (M, YO, BSQ, RHORAT , YN , IER  ) 

-  uses  a  Newton -R a phson  iteration  procedure  to 
calculate  YN ,  which  gives  the  M-th  normal  mode  eigenvalue 
I  see  Eqs.  29-311.  YO  is  an  initial  guess,  BSQ  and  RHORAT 
some  needed  input  variables,  and  IER  an  output  error 
parameter  in  case  the  method  fails  to  converge  to  the 
required  tolerance. 


3.  subroutine  ATTENU 

-  calculates  the  imaginary  part  of  the  wave  number 
for  the  modes  due  to  the  bottom  attenuation  using  a 
perturbation  formula.  Parameters  are  passed  using  labelled 
common . 
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subroutine  MODFUN  (NM,IO) 

-  calculates  and  prints  (on  logical  unit  number  10) 
the  NM  eigenfunctions  evaluated  at  the  12  depths  equally 
spaced  at  intervals  of  one  tenth  the  thickness  of  layer  one. 
This  subroutine  should  be  generalized  in  future  versions. 


5.  subroutine  WFUNCT  (N , NPT , DZ , 2 VA LS , WVALS , ICH  ) 

-  evaluates  and  stores  in  array  WVALS  the  mode 
function  for  the  N-th  mode  at  the  NPT  depths  stored  in  array 
ZVALS.  If  ICH  is  non-zero  the  array  ZVALS  is  calculated: 
ZVALS (I ) = (I -1  )  *DZ  .  The  mode  parameters  are  obtained  via 
labelled  common. 


6.  function  UN  ( Z , N , H , C 1 , C2 , RHORAT , CM  EGA , NM , AKM ,  ANM  ) 

-  evaluates  the  N-th  mode  function  at  depth  Z.  The 
other  parameters  are  required  for  calculation  of  the  mode 
function. 


7.  subroutine  PLOSS  (10) 

-  In  put/output  subroutine  for  propagation  loss 
calculations.  TO  is  the  logical  number  of  the  file  on  which 
the  results  are  written. 


8.  subroutine  EXCITE  (ZS,EX,NM) 

-  calculates  the  mode  excitations  (array  EX)  for 
source  depth  ZS  for  the  first  NM  modes. 


9.  function  PROPLO  ( R , ZR , ZS , COH I , RCOH I  ) 

-  for  range  R ,  receiver  depth  ZR,  and  source  depth 
ZS,  calculates  the  coherent  propagation  loss  PROPLO,  the 
incoherent  propagation  loss  COHI,  and  incoherent  propagation 
loss  RCOH I  with  cylindrical  and  near-field  spherical 
spreading  removed.  All  propagation  losses  are  in  dB 
referenced  to  unit  distance. 


10.  complex  function  VPHI  (R,ZR,ZS,VPHI2) 

-  evaluates  the  velocity  potential  VPHI  and  the 
square  of  the  incoherent  velocity  potential  VPHI2,  for  range 
R,  receiver  depth  ZR  and  source  depth  ZS. 


11.  complex  function  HANKC1  (RK) 


-  evaluates  the  asymptotic  Hankel  function  of  the 
zeroth  order  and  first  kind  for  large  values  of  RK. 


The  following  DEC-system  routines  are  used  [10]. 
They  are  not  essential  to  the  calculations  and  are  used  to 
uniquely  label  the  results. 


1.  subroutine  DATE  (DDMMYY) 

-  places  the  current  date  as  lef t- justi f ied  ASCII 
characters  into  a  two-word  array  DDMMYY.  The  date  is  of  the 
form  dd-mmm-yy,  where  dd  is  a  two-digit  day  (if  the  first 
digit  is  zero,  it  is  converted  to  a  blank),  mmm  is  a 
three-character  month  and  yy  is  a  two-digit  year;  e,g. 

3 1 -Mar-78 . 


2.  subroutine  TIME  (HHMM) 

-  returns  the  current  time  in  HHMM  as  left- justi f ied 
ASCII  characters  in  the  form  hh:mm,  where  hh  is  the  hours  in 
?iJ-hour  time  and  mm  is  the  minutes. 


A  .  7  A  o  c  u  r  a  c  y 

The  program  and  subroutines  are  written  in  single 
precision,  so  the  accuracy  is  machine  dependent  to  some 
extent.  The  DEC-20  has  36-bit  words,  or  about  8  digits 
accuracy  for  floating  point  numbers. 

With  the  tolerance  on  the  eigenvalue  set  at  1.0E-8, 
the  wave  numbers  seem  to  be  correct  to  8  significant 
figures,  the  mode  functions  seem  correct  to  6  significant 
figures,  and  for  the  non-attenuated  case  the  propagation 
loss  seems  correct  to  about  4  significant  digits.  More 
extensive  comparisons  are  made  in  section  3-3-  In  all  cases 


tried  so  far,  convergence  has  occurred  in  fewer  than  ?0 


iterations. 

Since  only  three  terms  are  used  in  the  expression 
for  the  Hankel  function,  propagation  loss  may  be  inaccurate 
at  short  ranges,  on  the  order  of  several  water  depths. 
However,  the  normal  mode  model  is  inaccurate  at  these  ranges 
anyway . 

A.  8  Memory  Requirements 

The  .EXE  file  requires  13  pages  or  6656  36-bit  words 
on  the  DEC-20. 

The  present  version  of  the  program  has  461  lines  of 
code,  a  total  of  12945  ASCII  characters. 

A. 9  Execution  times 

The  execution  time  depends  more  or  less  linearly  on 
the  number  of  modes  calculated.  The  number  of  modes  (see  Eq . 

26)  in  turn  depends  linearly  on  the  frequency  and  on  the 
thickness  of  layer  one.  The  number  of  modes  also  increases 
with  the  sound  speed  ratio  CRAT. 

For  the  calculation  of  Appendices  B  and  C,  with 
seven  modes  and  all  options  in  effect,  one  and  one-half  to 
three  CPU  seconds  execution  time  is  used  on  the  DEC-20.  The 
exact  execution  time  depends  on  the  system  load. 
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APPENDIX  B 


SAMPLE  TERMINAL  INPUT/OUTPUT 

All  user  entries  are  underlined.  The  symbol  <ret> 
indicates  that  the  user  must  press  the  "return"  key  to 
terminate  the  data  entry.  The  line  printer  output 
associated  with  this  example  is  given  in  Appendix  C. 
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IrtiECEDIM  PAOSl  SUM -HOT  F ILMD 


^EXECUTE  N0RM2L  <r£t> 
LINK:  Loading 

[ LNKXCT  N0RM2L  Execution] 

i4-Jun-79  09:30 


NORMAL  MODE  TWO  LAYER  MODEL 


WATER  DEPTH  ? 

TOO  <ret> 

FIRST  LAYER:  Cl ,RH01  ? 
i 500 .  1 ■  <r£t> 

SECOND  LAYER  RATIOS:  CRAT , RHORAT  ? 

2^  <r£t> 

C2.RH02  :  1800.000  2.0 

ATTENUATION  IN  DB/UNIT  LENGTH  //  1  KHZ  ? 

0.2  <r£t> 

FREQUENCY  ? 
iQO  <r£t> 

i 4-Jun-79  09:30 


NORMAL  MODE  TWO  LAYER  MODEL 


WATER  DEPTH  100.0 
Cl,  RH01  1500.00  1.0000 

C2 ,  RH02  1800.00  2.0000 

CRAT, RHORAT  1.2000  2.0000 

FREQUENCY  100.00 

ALPRAT  0.2000  DB/UNIT  DISTANCE  //  1  KHZ 

ATTN  0.00230  UNITS  OF  INVERSE  DISTANCE 

7  MODES  CALCULATED 
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MODE  KN  NORM  VPHASE  IM(KN) 

1  0.41787760  0.135884  1503.5947  0.23081E-05 

2  0.41483675  0.136394  1514.6163  0.85045E-05 

3  0.40965783  0.137016  1533.7642  0.17149E-04 

4  0.40220081  0.137577  1562.2010  0.27377E-04 

5  0.39229571  0.137972  1601.6451  0.3971  IE-04 

o  0.37974394  O.136IO1  1654.5847  0.57252E-04 

7  0.36434072  0.137507  1724.5356  0.94016E-04 


DO  YOU  WANT  A  COPY  OF  THIS  SUMMARY  ON  THE  LINE  PRINTER  ? 
TYPE  YES  OR  NO. 

Y  <rgt> 

CALCULATE  MODE  FUNCTIONS  ?  (Y/N). 

Y  <rgt> 

CALCULATE  PROPAGATION  LOSS  ?  (Y/N). 

Y  <rgt> 

CHANGE  DEFAULT  VALUES  ?  (Y/N). 

N  <rg.t> 

CALCULATE  EXCITATIONS  OF  EACH  MODE  ?  (Y/N) 

I  <r£ t> 

MORE  PROPAGATION  LOSS  CALCULATIONS  ?  (Y/N). 

Y  <r£t> 

CHANGE  DEFAULT  VALUES  ?  (Y/N). 

Y  <rgt> 

SOURCE  DEPTH,  RECEIVER  DEPTH  ? 

25 <r£t> 

RMIN , DELR , NRNG  ? 

. . .  ^r£t.> 

CALCULATE  EXCITATIONS  OF  EACH  MODE  ?  (Y/N) 
i  ^r£t> 

MORE  PROPAGATION  LOSS  CALCULATIONS  ?  (Y/N). 

N  <r£t> 

RUN  ANOTHER  DATA  SET  ?  (Y/N). 

JS.  <rjgt> 

END  OF  EXECUTION 

CPU  TIME:  2.09  ELAPSED  TIME:  1:13. 20 
EXIT 

e 


APPENDIX  C 

SAMPLE  LINE  PRINTER  OUTPUT 

The  output  shown  here  was  generated  by  the  terminal 
session  illustrated  in  Appendix  B. 
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l«eceoing  eaos  bunk -NOT  fi  la® 


14-J jn-79  09:30 


NORMAL  MODE  TWO  LAYER  MODEL 


WATER  DEPTH 
Ci,  RHOi 
CP,  RH02 
CRAT , RHORAT 


10C.0 

1500.00  1.0000 

1600.00  2.0000 

1.2000  2.0000 


FREQUENCY  100.00 


ALPRAT 

ATTN 


0.2000  DB/UNIT  DISTANCE  //  i  KHZ 
0.00230  UNITS  OF  INVERSE  DISTANCE 


7  MODES  CALCULATED 


MODE  KN 

1  0.41787760 

2  0.41483675 

3  0.40965783 

4  0.40220081 

5  0.39229571 

6  0.37974394 

7  0.36434072 


NORM 

0.135884 
0.136394 
0 . 1 37016 
0.137577 
0.137972 
0.138101 
0.137587 


VPHASE 

1503.5947 

1514.6163 

1533.7642 

1562.2010 

1601.6451 

1654.5847 

1724.5356 


IM(KN) 
0.23081E-05 
0.85045E-05 
0. 17 149E-04 
0.27377E-04 
O.397HE-04 
0.57252E-04 
0.94016E-04 


MODE  WAVE  FUNCTIONS 

DEPTHS 

10.00000 
50.00000  , 

90.00000  , 

MODES 

1 

0. 3878757E-01 , 
0.1348502 
0.6943669E-01 , 

O 

20.00000  , 

60.00000  , 

100.0000 

30.00000 
70.00000  , 
110.0000 

40.00000 
80.00000  , 
120.0000 

0.7434761E-01 , 

0.1340154 

0.3320670E-01, 

0.1037211 
0.122029 2  , 

0.1669175E-02, 

0.1244640  , 
0.9986882E-01 , 
0. 1678062E-03, 

C 

0.7480736E-01 , 
0.3228419E-01 , 

-O.H88813 

0.1251040 
-0.4568638E-01, 
-0.6273356E-01 , 

0.1544102 
-0.1086876  , 
-0.3334304E-02, 

0.9967681E-01 , 
-0.1360775 
-0.3544382E-03, 

J 

0.1050843  , 

-0.1290802 
0.1370046 

0.1348664 
-0.1180745 
0.8656157E-01 , 

0.6800471E-01 , 
-0.2245802E-01 , 
0.5072421E-02, 

-0.4758839E-01 , 
0.8925168E-01 , 
0.5943403E-03, 
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4 

0.1266855  ,  i 

-0.5761266E-01 ,  1 

-0.1230552  , 

5 

0.1372500 
0.1203002 
0.8347185E-01 , 

6 

0.1354279  ,  - 

0.7630017E-01 ,  - 
-0.2784262E-01 ,  - 
7 

0.1210085  ,  - 

-0.1085481  ,  - 

-0.3385370E-01 , 


PROPAGATION  LOSS 


SOURCE  DEPTH 

RANGE 

10000.00 

20000.00 

30000.00 

40000.00 

50000.00 

60000.00 

70000.00 

80000.00 

90000.00 

100000.0 


3.9860192E-01 , 
0.9257621E-01 , 
0.1046353 

0.2804273E-01 , 
0.7949435E-01 , 
0.1178102  , 

0.5303766E-01 , 
0.1278222 
0.1271950 

0.1151735 
0.2269905E-01 , 
0.1333989 


-0.4962995E-01 , 
0.1298128 
-0.7094869E-02, 

-0.1315203 

-0.1040581 

0.9832619E-02, 

-0.1146567 

-0.26241 15E-01 , 
-0. 1425772E-01 , 

-0-11 3886 1 E-0 1 , 
0.1301526 
0.2348356E-01 , 


-0.1375083 

0.8664670E-02, 

-0.9621448E-03, 

-0.5491478E-01 
-0.1007553 
0. 1641291E-02 

0.9794070E-01 
0.1 380990 
-0.3196392E-02 

0.1260130 

-0.1011776 

0.8268099E-02 


MODE  EXCITATIONS  SOURCE  DEPTH  =  50. 

1  0. 13485018 

2  0. 32284 186E-01 

3  -0.12908018 

4  -0.57612656E-01 

5  0.12030022 

6  0.76300173E-Q1 

7  -0.10854806 


50.00 

RECEIVER  DEPTH 

50.00 

P.L. 

(DB  RE  UNIT  DISTANCE) 

COH 

INCOH 

INCOH-GEOM 

58.37791 

,  60.02504 

0.2504134E-01 , 

64.97957 

,  64.23382 

1.223522 

65.42920 

,  66.82623 

2.055018 

70.69642 

,  68.70615 

2.685551 

76.35365 

,  70.17996 

3.190257  , 

75.91258 

,  71.39203 

3.610515 

74.76040 

,  72.42286 

3.971880 

74.88633 

,  73.32203 

4.291131 

74.89706 

,  74.12217 

4.579747 

74.40432 

,  74.84582 

4.845824 
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PROPAGATION  LOSS 


MODE  EXCITATIONS 

1  0.89975154E-01 

2  0.13542200 

3  0.11195862 

4  0.29491923E-01 

5  -0.69696217E-01 

6  -0.13222866 

7  -0.12361683 


SOURCE  DEPTH  =  25.00 


SOURCE  DEPTH 

25.00 

RECEIVER  DEPTH 

50.00 

RANGE 

P.L. 

(DB  RE  UNIT  DISTANCE) 

C0H 

INCOH 

INCOH-GEOM 

10000.00 

> 

67.96269 

,  62.09322 

2.093222 

20000.00 

1 

66.87508 

,  66.65010  , 

3.639804 

30000.00 

7 

71.06421 

,  69.42616 

4.654949 

40000.00 

t 

73.10997 

,  71.44136 

5.420760 

50000.00 

7 

69.76201 

,  73.03228  , 

6.042584 

60000.00 

7 

70.69862 

,  74.34834 

6.566830 

70000.00 

7 

72.88694 

,  75.47003 

7.019045 

80000.00 

7 

75.62021 

,  76.44690 

7.416999 

90000.00 

7 

78.87690 

,  77.31236 

7.769934 

100000.0 

7 

80.43400 

,  78.09026 

8.090258 
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APPENDIX  D 


PROGRAM  LISTING 


PROGRAM  NORM2L 
C 

C  TWO  LAYER  NORMAL  MODE  PROGRAM 

C 

C  AUTHOR:  DALE  D.  ELLIS 

C  DEFENCE  RESEARCH  ESTABLISHMENT  ATLANTIC 

C  DARTMOUTH,  NOVA  SCOTIA,  CANADA 

C 

C  DATE  WRITTEN:  APRIL  1978 

C  LAST  MODIFICATION:  JUNE  1978 

C  (  SOME  COSMETIC  CHANGES  SINCE  THEN) 

C  VERSION  :  5.0  MOD  41 

C  COMPUTER:  DEC-20/40 

C  OPERATING  SYSTEM:  TOPS-20 

C 

C  LANGUAGE:  FORTRAN  IV,  WITH  A  FEW  DEC-SYSTEM  EXTENSIONS 

C 

c 

c 

COMMON  /ENVIRN/  H ,C1 , RH01 ,C2 , RHORAT ,ATTN .OMEGA ,AK1 ,AK2 
COMMON  /EIGENF/  NM,AKM(50) , ANM(50) ,ATN(50 ) 

COMMON  /CNTLNM/  MAX ,TOL ,NMIN ,NMAX , VMIN , VMAX 
COMMON  /CNTLIO/  II, 10, IT 
DIMENSION  DDMMYY (2 ) 

DATA  H , C 1 , RHO 1 , CHAT , RHORAT , ALPRAT , FREQ/ 100 . , 1500 . , l . , 1 . 2 

i,  2.0,0.2,100.  / 

DATA  II, 10, IT  /5 , 3 ,5/ 

DATA  TOL  /1.E-8/,  MAX/20/,  NMIN/1/,  NMAX/50/ 

DATA  NCALL/O/ 

C 

PI2  =  8.»ATAN(1.) 

ISAVE=IO 

DBNP=ALOG(10. )/20. 

CALL  DATE  (DDMMYY) 

CALL  TIME  (HHMM) 

WRITE  (IT, 500)  DDMMYY, HHMM 
500  FORMAT  ( 1X,2A5,5X,A5) 

C 

C  INPUT  DATA 

C 


2  CONTINUE 
C 


WRITE 

(IT, 400) 

400 

r * 

FORMAT 

( /// , 1 1 X , 28H  NORMAL  MODE  TWO  LAYER  MODEL  ) 

WRITE 

(IT, 401) 

401 

FORMAT 

( /// ,  14H  WATER  DEPTH  ?) 

r> 

READ  (II,*)  H 

L 

4  WHITE 

(IT, 402) 

402 

FORMAT 

(24H  FIRST  LAYER:  C1.RH01  7  ) 

READ  (11,*)  C 1 , RHO 1 

IF ( ABS ( RHO 1 - 1 . )  .LT.0.03)  GO  TO  5 
WRITE(IT,412) 

412  FORMAT  (•  ARE  YOU  SURE  ABOUT  RHOI  ?  (Y/N).') 

READ  (11,302)  YESNO 

IF(YESNO.NE. 1HY)  GO  TO  4 
5  WRITE  (IT, 403) 

403  FORMAT  (36H  SECOND  LAYER  RATIOS:  CRAT.RHORAT  ?  ) 

READ  (II,*)  CRAT.RHORAT 

IF  (CRAT.GT. 1 . )  GO  TO  6 
WRITE  (IT, 417) 

417  FORMAT  ('  CRAT=C2/C1  MUST  BE  GREATER  THAN  1.0.  ENTER  NEW  VALUE.’) 
GO  TO  5 

b  RH02  =  RH01*RH0RAT 
C2  =  CRAT  *  Cl 
WRITE  (IT,  404)  C2.RH02 

404  FORMAT  (22X , 10HC2 ,RH02  :  ,F10 . 3 ,F10 . 1 , /) 

C 

WRITE  (IT, 406) 

406  FORMAT  ('  ATTENUATION  IN  DB/UN1T  LENGTH  //  i  KHZ  ?’) 

READ  (II,*)  ALPRAT 
C 

WRITE  (IT, 405) 

405  FORMAT  (12H  FREQUENCY  ?) 

READ  (II,*)  FREQ 

OMEGA  =PI2*FREQ 
IF  (FREQ.LT.O. )  WRITE  (IT, 413) 

413  FORMAT  ('  NEGATIVE  FREQUENCY.  EXECUTION  SUPPRESSED.’) 

IF  (FREQ.LT.O.)  GO  TO  25 

C 

HOMC=  H*0MEGA/C1 
AK1  =  OMEGA/C  1 
AK2  =  0MEGA/C2 

ATTN=DBNP*(ALPRAT*FREQ/1000. ) 

C 

CALL  NMODES  ( H , HOMC , CRAT , RHO 1 , RHOR AT , NM , AKM , ANM ) 

C 

IF  (NM.LE.NMAX)  GO  TO  8 
WRITE  (IT, 419)  NM,NMAX 

419  FORMAT (//,'  ***  WARNING  ***  THERE  ARE ’, 15 , '  MODES',/ 

1,'  ONLY ',15,'  MODES  HAVE  BEEN  CALCULATED.’) 

NM=NMAX 

8  IF  (NM.GT.O. )  CALL  ATTENU 
10=IT 

1C  =  0 

9  CONTINUE 

IF(NCALL.NE.O)  WRITE  (10,499) 

499  FORMAT(IHl) 

WRITE  (10,500)  DDMMYY ,HHMM 
WRITE  (10,400) 

WRITE  (10,501)  H, Ci, RHO’, C2.RH02, CRAT.RHORAT, FREQ 
1,  ALPRAT, ATTN 
WRITE  (10,507)  NM 
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FORMAT  (/,I5, '  MODES  CALCULATED') 

IF  (NM.LT.i)  GO  TO  13 
WRITE  (10,510) 

FORMAT  (//,6H  MODE , 7X , 2HKN , 1 3X , 4HNORM, 10X , 6HVPHASE ,9X , 

1  6HIM(KN)  ) 

DO  11  1=1, NM 
VPHASE=OMEGA/AKM(I ) 

WRITE  (10,301)  I , AKM(I ) , ANM(I ) , VPHASE , ATN(I ) 

11  CONTINUE 

FORMAT  (I5,G17.8,G15.6,G16.8,G13.5) 

FORMAT  (12H-WATER  DEPTH, F10. 1 ,/, 12H  Cl,  RH01  , F1 1.2, 

1  F1 1 . 4 , / , 12H  C2 ,  RH02  ,F1 1 .2 ,F1 1 .4 , /12H  CRAT.RHORAT, 

2  F13.4,F9.4,//,12H  FREQUENCY  ,F11.2  , 

3  //,7H  ALPRAT.F18.4, '  DB/UNIT  DISTANCE  //  1  KHZ’, 

4  /,5H  ATTN, 1 1X, FI 0.5,'  UNITS  OF  INVERSE  DISTANCE') 

13  IF  (IC.EQ.1)  GO  TO  15 

WRITE  (IT, 407) 

FORMAT  (//,'  DO  YOU  WANT  A  COPY  OF  THIS  SUMMARY  ON  THE  LINE 
1'  PRINTER  ?' ,/, 10X, 'TYPE  YES  OR  NO.') 

READ  (IT, 302)  YESNO 
FORMAT  ( A1  ) 

IF  (YESN0.NE.1HY)  GO  TO  15 

IO=ISAVE 

IC=  1 

NCALL=NCALL+1 
GO  TO  9 

15  CONTINUE 

IF  (NM.LT.I)  GO  TO  25 
IO=ISAVE 
WRITE  (IT, 408) 

FORMAT  ('  CALCULATE  MODE  FUNCTIONS  ?  (Y/N).') 

READ  (IT,  302)  YESNO 
IF  (YESNO. NE. 1HY)  GO  TO  20 
CALI.  MODFUN  (NM,IO) 

20  WRITE  (IT, 409) 

FORMAT (  '  CALCULATE  PROPAGATION  LOSS  ?  (Y/N).') 

READ  (IT, 302)  YESNO 
IF  (YESNO. NE. 1HY  )  GO  TO  25 
22  CALL  PLOSS  (10) 

WRITE  (IT, 4 15) 

FORMAT  ('  MORE  PROPAGATION  LOSS  CALCULATIONS  ?  (Y/N).') 

READ  (IT, 302)  YESNO 
IF  ( YESNO. EQ.1HY)  GO  TO  22 
25  WRITE  (IT, 410) 

1  FORMAT ('  RUN  ANOTHER  DATA  SET  ?  (Y/N).’) 

READ  (IT, 302)  YESNO 
CALL  TIME  (HHMM) 

IF  (YESNO. EQ.1HY)  GO  TO  2 
END 


o  o 


SUBROUTINE  NMODES  ( H , HOMO , CRAT , RHO 1 , RHORAT , NM , AKM , ANM ) 

C 

NORMAL  MODE  CALCULATION  FOR  TWO  LAYERS 

C  WRITTEN  BY  DALE  D  ELLIS  APRIL  1978 

C 

COMMON  /CNTLNM/MAX , TOL , NMIN , NMAX , VMIN , VMAX 
COMMON  /CNTLIO/  II, 10, IT 
LOGICAL  EXTRA 

DIMENSION  AKM ( NM ) , ANM ( NM ) 

DATA  PI  /3. 14159265  / 

DATA  EXTRA  /.FALSE./ 

BSQ  =  HOMC  **2  *  (1.-1 ./CRAT**2  ) 

B  =  SQRT  (BSQ) 

NM  =  B/PI  +  0.5 
IF  (NM.LT.1)  GO  TO  40 
NMX  =  MINO  (NM.NMAX ) 

YO  =  (NMIN-i .25 )*PI 
DO  20  M  =NMIN , NMX 
YO  =  YO+PI 

IF (EXTRA )  WRITE (IT ,410 )  M 
410  FORMAT  (5HOMODE.I5) 

IM  =  M 

CALL  MODE  ( IM.YO, BSQ, RHORAT, YN.IER) 

RKNSQ=  (H0MC**2-YN**2 ) /H**2 
RKN  =  SQRT  (RKNSQ) 

AKM(M) =RKN 
GM1H  =  YN 

GM2H  =  SQRT  (  BSQ-YN**2) 

AM2  =  2./H/RH01/(1.-SIN(2.*GMlH)/(2.*GMiH)  +(SIN  (GM1H) ) 
1 **2 / ( RH0RAT*GM2H ) ) 

ANM(M) =SQRT ( AM2 ) 

IF  (EXTRA)  WRITE  (IT,*)  YN , AKM(M) ,ANM(M) 

20  CONTINUE 

C 

RETURN 

C 

40  WRITE (IT, 420) 

420  FORMAT  PYH  NO  TRAPPED  MODES  ) 

RETURN 

END 
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SUBROUTINE  MODE  (M, YO , BSQ , RHORAT , YN , IER ) 


C 

C  OBTAINS  YM  FOR  M-TH  MODE  USING  NEWTON-RAPHSON  PROCEDURE 

C 

C  WRITTEN  BY  DALE  D  ELLIS  APRIL  1978 

C 

LOGICAL  EXTRA 

COMMON  /CNTLNM/  MAX , TOL , NMIN , NMAX , VMIN , VMAX 
COMMON  /CNTLIO/  II ,10, IT 
DATA  PI  /3. 14159265358979/ 

DATA  EXTRA  /.FALSE./ 

IER:0 

B=SQRT(BSQ) 

R21B= (RH0RAT**2-1 . )/BSQ 
PIM=PI*M 

YMAX=AMIN1 (PIM,B) 

YMIN=PI*(M-0 . 5 ) 

YM=Y0 

IF  ( YM. GT . YMAX . OR . YM. LT . YMIN )  YM= (YMIN+YMAX) /2 . 

C 

DO  10  1=1, MAX 

SRBY  :  SORT  (BSQ-YM*»2) 

FY  =  YM-PIM  +  ATAN  (RHORAT*YM/SRBY ) 

DF  =  1.  +  RHORAT/SRBY/0.  +  R21B*Y»»2) 

YN  =  YM-FY/DF 
DIF=  YN-YM 

IF  (EXTRA)  WRITEUT,*)  YM,FY,DF,YN,DIF 
IF  (ABS(DIF/YN).LT.  TOL)  GO  TO  15 
YM=YN 

IF  (YM.LT.YMIN )  YM=YMIN 
IF  (YM.GT.YMAX)  YM=YMAX 
10  CONTINUE 

C 

WRITE  (IT, 410)  M, MAX, DIF 

410  FORMAT  (34H  CONVERGENCE  NOT  ACHIEVED  FOR  MODE  ,I3,/,10X, 
1  5HAFTER , 13 , 19H  ITERATIONS  DIF  =,  E12.4) 

IER:  1 
RETURN 

C 

15  CONTINUE 
HETURN 
END 


SUBROUTINE  ATTENU 


C  BOTTOM  ATTENUATION  USING  INGENTO'S  PERTURBATION  FORMULA 
C 

C  WRITTEN  BY  DALE  D  ELLIS  APRIL  1978 

C 

COMMON  /ENVIRN/  H ,C1 ,RH01 ,02 .RHORAT , ATTN .OMEGA, AK1 ,AK2 
COMMON  /EIGENF/  NM,AKM(50) ,ANM(50) ,ATN(50) 

C 

IF  (ATTN.EQ.O. )  GO  TO  20 
C 

CONST =0 . 5*ATTN#0MEGA*RH0  VC2/RH0RAT 

C 

DO  10  1=1, NM 

ATN(l)  =  CONST*( ANM( I ) *S1N (H*SQRT ( AK1 **2-AKM(I )**2 ) ) )**2 
1  /  (AKM(I )*SQRT(AKM(I )*#2-AK2**2)  ) 

10  CONTINUE 
RETURN 
C 

C  NO  ATTENUATION 

20  DO  21  1=1, NM 

21  ATN (I ) =0 . 

RETURN 

END 
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SUBROUTINE  MODFUN  (NM,IO) 

INPUT/OUTPUT  ROUTINE  FOR  MODE  FUNCTION  CALCULATIONS 

WRITTEN  BY  DALE  D  ELLIS  APRIL  1978 

C 

COMMON  /ENVIRN/  H,C1 ,RH01 ,C2,RHORAT, ATTN, OMEGA, AK1 ,AK2 
DIMENSION  ZVALS(IOO) ,WVALS(100) 

WRITE  (10,430) 

430  FORMAT (////, 1 0X, 'MODE  WAVE  FUNCTIONS',/) 

NPT=12 

IN= 1 

DZ=H/10 . 

CALL  WFUNCT  ( 1 , NPT ,DZ , ZVALS ,WVALS ,0  ) 

WRITE  (10,431) 

431  FORMAT  ('  DEPTHS') 

WRITE  (10,*)  (ZVALS (I ) ,1  =  1 ,NPT ) 

WRITE  (10,432) 

432  FORMAT  (/,'  MODES') 

WRITE  (10,202)  IN 

202  FORMAT  (lj) 

WRITE  (10,*)  (WVALS(l) ,1=1 ,NPT) 

IF  (NM.LT.2)  RETURN 
DO  10  IN=2,NM 
INN=IN 

CALL  WFUNCT  (INN, NPT, DZ, ZVALS, WVALS.O) 

WRITE  (10,202)  INN 

WRITE  (10,*)  (WVALS(I ) , 1= 1 ,NPT ) 

10  CONTINUE 
RETURN 
END 
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SUBROUTINE  WFUNCT  (N , NPT , DZ , ZVALS ,WVALS , ICH ) 


EVALUATES  MODE  FUNCTIONS  AT  A  NUMBER  OF  DEPTHS 

WRITTEN  BY  DALE  D  ELLIS  APRIL  1978 

COMMON  /ENVIRN/  H ,C1 , RHOI , C2 , RHORAT , ATTN .OMEGA ,AK1 ,AK2 
COMMON  /EIGENF/  NM,AKM(5Q) ,ANM(50) ,ATN(50) 

DIMENSION  ZVALS(NPT) .WVALS(NPT) 

IF  (ICH.NE.O)  GO  TO  10 
DO  4  1=1 , NPT 
4  ZVALS  (1 ) =DZ*I 
10  CONTINUE 

DO  14  1=1, NPT 

WVALS(I)  =  UN  ( ZVALS (I),N,H,C1,C2, RHORAT , OMEGA , NM , AKM , ANM  ) 
14  CONTINUE 
RETURN 
END 


FUNCTION  UN  (Z, N.H, Cl .C2.RH0RAT, OMEGA, NM.AKM.ANM  ) 

N-TH  NORMAL  MODE  EVALUATED  AT  DEPTH  Z 

WRITTEN  BY  DALE  D  ELLIS  APRIL  1978 

DIMENSION  AKM(NM) ,ANM(NM) 

IF  (N.GT.NM)  GO  TO  13 

GM1=SQRT  (  (OMEGA/C1  )*»2  -  AKM(N)*«2  ) 

IF  (Z.GT.H)  GO  TO  4 
UN  =  ANM(N )*SIN  (GM1*Z) 

RETURN 

4  UN= 1 . /RHORAT *ANM ( N ) *EXP ( - ( Z-H ) *SQRT ( AKM ( N ) **2- ( OMEGA/C2 ) «*2 ) ) 
1  *  SIN  (GM1*H) 

RETURN 
13  UN=0. 

RETURN 

END 


SUBROUTINE  PLOSS  (10) 


C 

C  INPUT/OUTPUT  ROUTINE  FOR  PROPAGATION  LOSS  CALCULATIONS 

C 

C  WRITTEN  BY  DALE  D  ELLIS  APRIL  1976 

C 

COMMON  /ENV1RN/  H ,C i , RH01 ,C2 .RHORAT , ATTN .OMEGA ,AKi , AK2 
COMMON  /EIGENF /  NM , AKM ( 50 ) ,  ANM  ( 50 ) ,  Al'N  ( 50 ) 

COMMON  /CNTL10/  II, ION, IT 
DIMENSION  EX (50) 

WRITE  (10,440) 

440  FORMAT  (///,'  PROPAGATION  LOSS'/) 

ZS=H/2 . 

ZR=H/2 . 

RMINr 10000. 

DELR= 10000. 

NRNG=10 

WDEPTH=  iO.*ALOGiO(H) 

C 

WRITE  (IT, 434) 

434  FORMAT  ('  CHANGE  DEFAULT  VALUES  ?  (Y/N).') 

READ(1T,435)  YESNO 

435  FORMAT  (A1) 

IF ( YESNO . NE . 1 HY )  GO  TO  6 
WRITE  (IT, 436) 

436  FORMAT  ('  SOURCE  DEPTH,  RECEIVER  DEPTH  ?') 

READ  (IT,*)  ZS , ZR 

WRITE  (IT, 437) 

437  FORMAT  ('  RMIN , DELR , NRNG  ?') 

READ  (IT,*)  RMIN, DELR, NRNG 

C 

6  CONTINUE 

WRITE  (IT, 436) 

438  FORMAT  ('  CALCULATE  EXCITATIONS  OF  EACH  MODE  ?  (Y/N)') 

READ  (IT, 435)  YESNO 

IF  (YESNO.NE.iHY)  GO  TO  8 
CALL  EXCITE  (ZS,EX,NM) 

WRITE  (10,439)  ZS 

439  FORMAT  (20X,'MODE  EXCITATIONS' , 10X, 'SOURCE  DEPTH  =',G11.4) 
WRITE  (10,444)  (I ,EX (I ) , 1=1 ,NM) 

444  FORMAT  ( 15X ,  15, d7 .8) 

6  CONTINUE 
C 

WRITE  (10,441)  ZS,ZR 

441  FORMAT  ( / , 1 3H  SOURCE  DEPTH ,F10.2, 10X, 14HRECEIVER  DEPTH, 

1  F10 .2 , // ,5X , 5H RANGE , 1 3X , ' P .L .  (DB  RE  UNIT  DISTANCE)’) 

WRITE  (Ju,442) 

442  FORMAT  (2  * X , ' COH ' , 12X , ' INCOH ' , 9X , ' INCOH-GEOM' ) 

DO  10  IRrl , NRNG 

R=RMIN+(IR-1 )*DELR 
PL=PROPLO  ( R , ZR , ZS , COHI , RCOHI ) 

GCOHI=RCOHI-WDEPTH 
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WRITE  (10,*)  R , PL , COHI ,GCOHI 
10  CONTINUE 
RETURN 
END 


SUBROUTINE  EXCITE  (ZS,EX,NMD) 

C 

C  CALCULATES  MODE  EXCITATIONS  (I.E.  MODE  FUNCTION  AT  SOURCE  DEPTH) 
C 

C  WRITTEN  BY  DALE  D  ELLIS  APRIL  1978 

C 

DIMENSION  EX (NMD) 

DO  10  I  = 1 , NMD 
II=I 

CALL  WFUNCT  (II , 1 ,DZ,ZS,EX(II) ,  1 ) 

10  CONTINUE 
RETURN 
END 


FUNCTION  PROPLO  ( R , ZR , ZS , COHI , RC0H1 ) 

C 

C  PROPAGATION  LOSS  IN  DB 

C 

C  WRITTEN  BY  DALE  D  ELLIS  APRIL  1978 

REVISED  JUNE  1979  -  QUANTITIES  MADE  POSITIVE 

COMPLEX  VPHI.V 
V  =  VPHI  ( R , ZR , ZS , VCOHI ) 

VMOD  =  V*CONJG  (V) 

PROPLO=  -10.*AL0G10  (VMOD) 

COHI  =  -10 , *ALOG 10 (VCOHI ) 

RCOHI  =  -10,*ALOG10(R*VCOHI) 

RETURN 
END 
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COMPLEX  FUNCTION  VPHI  (R ,ZR ,ZS,COHI ) 

C 

C  VELOCITY  POTENTIAL  FOR  PROPAGATION  LOSS  CALCULATIONS 
C  ALLOWS  RECEIVER  DEPTH  ZR  TO  BE  GREATER  THAN  WATER  DEPTH  H 

C 

WRITTEN  BY  DALE  D  ELLIS  APRIL  1978 

COMPLEX  HANKOi ,HK,TERM 
DOUBLE  PRECISION  SUMS, SUMI 

COMMON  /ENVIRN/  H , C 1 , RHO 1 , C2 , RHORAT , ATTN , OMEGA , AK 1 , AK2 
COMMON  /EIGENF/  NM,AKM(50) ,ANM(50) ,ATN(50) 

DATA  PI/3. 141 59265 / 

C 

SUMR=0 . DO 
SUMI=0.D0 
COHI=0 . 

DO  10  1=1 ,NM 
HK=HANKO 1 (AKM(I)*R) 

GMl =SQRT  (AK1**2-AKM(I)**2) 

UNS  =SIN  (ZS*GM1)  *ANM(I) 

ALOSS=EXP (-ATN (I )*R) 

IF(ZR.GT.H)  GO  TO  6 
UNR  =SIN  (ZR«GM1)  *ANM(I) 

GO  TO  8 

6  GM2=  SQRT  (AKM(I)**2-AK2**2) 

UNR  =1./  RHORAT *ANM ( I ) *EXP ( - ( ZR-H ) *GM2 )  »  SIN  (GMl *H ) 

8  SUMR  =  SUMR  +  UNS«UNR*ALuSS»REAL(HK) 

SUMI  =  SUMI  +  UNS*UNR*ALOSS«AIMAG(HK) 

TERM= ( UNS*UNR*ALOSS ) *HK 
COHI=TERM*CONJG (TERM)  +  COHI 
10  CONTINUE 
C 

VPHI  =  CMPLX (0 . , PI )*CMPLX(SNGL(SUMR ) ,SNGL (SUMI ) ) 

C0HI=PI»*2  *  COHI 

RETURN 

END 
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COMPLEX  FUNCTION  HANKOKZ) 


C 

C  ASYMPTOTIC  HANKEL  FUNCTION  OF  FIRST  KIND  AND  ZEROTH  ORDER 
C 

C  WRITTEN  BY  DALE  D  ELLIS  APRIL  1978 

C 

DATA  PI  /3. 14159265/ 

CHI  =  Z  -  PI/4. 

P  =  1.  -  (9./128J/Z **2 
Q  =  -1./(8*Z) 

C 

HANKOI  =  SORT (2. /(PI *Z) )*CMPLX(P ,Q)*CMPLX(COS(CHI ) , SIN (CHI ) ) 

RETURN 

END 
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APPENDIX  E 


FUTURE  MODIFICATIONS  AND  EXTENSIONS 


Some  suggested  improvements  and  extensions  to  the 
present  program  are  listed  below. 


1.  Make  optional  instructions  available  to  the  user  at 
the  beginning  of  program  execution. 

2.  Type  the  current  value  of  the  data  item  to  be  entered. 
For  example,  the  first  question  might  read: 

'WATER  DEPTH  100.  ?' 

3.  Add  a  subroutine  CHOICE  to  allow  the  user  to  modify 
some  of  the  control  parameters  in  the  COMMON  block 
/CNTLNM/.  For  the  option  of  starting  the  calculations  at 
a  mode  number  greater  than  1,  some  subroutines  would  need 
modification . 

*4.  Improve  the  Hankel  function  to  maintain  accuracy  near 
the  source.  Type  a  warning  message  when  the  function  is 
not  sufficiently  accurate. 

5.  Generalize  MODFUN  to  handle  an  arbitrary  number  of 
specified  depths.  The  subroutines  called  by  MODFUN  are 
sufficiently  general  to  handle  this  possibility. 

6.  Add  the  possibility  of  the  optional  printout  appearing 
on  both  the  line  printer  and  terminal. 

*7.  Calculate  the  group  velocity,  cutoff  frequency,  and 
equivalent  ray  angle  for  each  mode. 

*8.  Optionally  display  the  mode  functions  and  propagation 
loss,  either  as  a  line  printer  plot  or  a  proper  graph. 


*  These  items  have  been  included  in  a  newer  version  of 
N0RM2L  -  PF.KRIS  -  which  is  still  undergoing  development. 
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It  is  hoped  to  try  some  new  ideas  using  N0RM2L  or  a 


revised  version  of  it. 


1.  Incorporate  range  dependent  normal  mode  theory,  at  least 
in  the  adiabatic  approximation. 

2.  Include  the  effects  of  scattering  from  a  rough  surface 
or  bottom. 

3.  Add  the  effects  of  continuous  modes.  This  is  necessary 
if  pressure  or  propagation  loss  is  needed  near  the  source. 

4.  Create  a  similar  program  in  which  the  bottom  layer  will 
also  support  shear  waves. 

5.  Already  parts  of  the  program  have  been  modified  and  used 
in  a  problem  involving  mode  enhancement  [4],  and  for 
debugging  a  general  normal  mode  program  [33-  These  uses 
were  not  foreseen  at  the  time  of  writing  the  program. 


APPENDIX  F 


USER  COMMENTS 

This  space  is  provided  for  the  user  to  add  his  notes 
or  comments.  Any  suggestions  for  improvements  or  additions 
will  be  considered  by  the  author  for  implementation  in 
future  versions  of  N0RM2L. 
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